astro-ph/9510102: Submitted to the Astrophysical Journal Letters 



Producing Mega-pixel CMB Maps from Differential Radiometer Data 
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ABSTRACT 

A major goal of cosmology is to obtain sensitive, high resolution maps of the 
Cosmic Microwave Background (CMB) anisotropy. Such maps, as would be produced 
by the recently proposed Microwave Anisotropy Probe {MAP), will contain a wealth 
of primary information about conditions in the early universe. To mitigate systematic 
effects when observing the microwave background, it is desirable for the raw data to 
be collected in differential form: as a set of temperature differences between points in 
the sky. However, the production of large (mega-pixel) maps from a set of temperature 
differences is a potentially severe computational challenge. We present a new technique 
for producing maps from differential radiometer data that has a computational cost 
that grows in the slowest possible way with increasing angular resolution and number 
of map pixels. The required central processor (CPU) time is proportional to the 
number of differential data points and the required random access memory (RAM) is 
proportional to the number of map pixels. We test our technique, and demonstrate its 
feasibility, by simulating one year of a space-borne anisotropy mission. 
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Introduction 



Experiments to measure the anisotropy of the Cosmic Microwave Background Radiation 
need to measure signals AT of 0(10 /xK) in the presence of the large isotropic background with 
To = 2.73 K and system noise that is usually 0(100 K). The difficulty inherent in stabilizing 
the sensitivity of an instrument to much better than 1 part per million has led all CMB 
anisotropy experiments to use differential radiometers, which chop rapidly between the spot whose 
temperature is to be measured and one or more reference sources. The reference sources are 
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usually other parts of the sky, although occasionally stabilized loads have been used. The resulting 

measurements of AT can be analyzed directly in terms of model predictions, but for high resolution 
experiments covering a large fraction of the sky, the number of independent data points can make 
analysis extremely cumbersome. It is far more desirable to reduce the data by producing a map 
of the microwave temperature. For a given angular resolution, a map provides the most complete 
form of data reduction possible without loss of information, which in turn permits a full range of 
statistical tests to be performed on the data. The problem of converting a large set of differential 
sky observations into a map becomes increasingly challenging as the angular resolution, number 
of map pixels, and data sampling rate all increase. The methods employed by COBE DMR can 
not be directly adapted to the high resolution case because of limitations imposed by computer 
RAM. We present a new technique for producing maps which is equivalent to the COBE sparse 
matrix technique, and which has a computational cost that grows in the slowest possible way with 
increasing angular resolution and number of map pixels: the required CPU time is proportional 
to the number of differential data points and the required RAM is proportional to the number of 
map pixels. We demonstrate the feasibility of this technique by simulating a one- year differential 
mission that produces sky maps with 1,572,864 square pixels of size 0?16. 



2. The DMR Experience 

The only full-sky CMB maps produced to date use the differential data from the COBE DMR. 
The DMR data analysis required the construction of sky maps using the 6 x 10^ differences that 
were collected each year in each of the six channels. Since there are only about 10^ beam areas on 
the sky, the system relating differences to sky temperature is highly over-determined. The COBE 
team chose to analyze the data using 6,144 pixels to cover the sky. These all have approximately 
the same area 9"^^^ and are arranged in a square grid of 32 x 32 pixels on each of the 6 faces of a 
cube. Within each pixel it was assumed that the temperature is constant, so the basic problem 
can be represented as a least-squares system with 24 x 10^ equations (after 4 years of data taking) 
in 6,144 variables. It is also desirable to extend this system to additionally account for systematic 
signals in the data due, for example, to the magnetic sensitivity of the ferrite switches used in 
the DMR instrument (Kogut et al. 1992). After calibration and baseline subtraction, one can 
represent the predicted output of the DMR instrument, in mK, as 

Sit) = vm^ (1) 

where t is the observation time. V{t) is a vector of length Npi^eis + 3 given by 

V{t) = [0, . . . , 0, +1, 0, . . . , 0, -1, 0, . . . , B,{t),By{t),B,{t)] (2) 

with the -1-1 in the pixel p{t) that contains the plus horn line-of-sight at time t, the -1 in the 
pixel m{t) that contains the minus horn line-of-sight, and Bi{t) being the ith component of the 
magnetic field in spacecraft-fixed coordinates. The parameter vector X is given by 

X = [To, ... , r6i43, M^, My, M,] (3) 
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where Tj is the temperature of the i*" pixel in mK, and and the Mj are the susceptibihties of the 

ferrite switches to external magnetic fields, in mK Gauss"^. For a constant noise per observation 
the least squares solution for the map is (Lineweaver et al. 1994; Wright et al. 1994a; Wright et 
al. 1994b) 

A = J2V{tfV{t) 

t 

t 

X = lim,_o+{A + eiy^ B (4) 

Note that A is sparse, symmetric, and singular. The main diagonal of A has elements An = N°'^^, 
the number of times each pixel was observed. The off-diagonal elements of A indicate the number 
of times a given pixel-pair ij was observed, either i = p{t) and j = m{t) or vice-versa. Since the 
angle between the DMR horns was fixed at 60°, only pixel pairs separated by 60° i 9pix could ever 
be observed. Thus, while A has 38 x 10^ elements, all but 1.8 x 10^ of these elements are zero, and 
since A^j = Aj^, only 9 x 10^ elements need to be kept. The limit e ^ 0"*" regulates the singularity 
of A, and gives a map with zero mean. 

The matrix A is sparse, but inverting A + el does not preserve this sparseness. Thus an 
iterative solution to the system {A + eI)X = B is found using algorithms that only require one to 
be able to find the product of A and a vector Y. If we let D be a diagonal matrix whose diagonal 
elements are AT?''*, then a very simple iterative scheme for computing X is 

Xin+l) ^ ^(n) ^ _ (5) 

Since diff^erential observations do not constrain the mean of the map, the mean should be 
adjusted to zero after each iteration. With the DMR scan pattern this simple iteration scheme 
converges quite rapidly. The computational cost of this algorithm for both RAM and CPU time is 
proportional to the number of non-zero off-diagonal elements in A. 



3. Scaling for Smaller Beams 

The number of off-diagonal elements can be calculated by considering the number of pixels in 
the reference ring of a given pixel. The radius of the reference ring is the chop angle, ^c/iop- The 
width of the reference ring is about twice the pixel size, 29pix- The solid angle of the the reference 
ring is then AnOpix sin Ochopi a-^d the number of pixels in the reference ring is Air sin 9 chop/ 9 pix- 
Since the total number of pixels is Npix = 47r/0pj^, the number of non-zero off-diagonal elements 
in the upper triangle of A is N^ff « Sn^ sin 9 chop/ f^^x- proposed new generation of 

anisotropy experiments, COBRAS/SAMBA, FIRE, PSI and MAP, the beam sizes will be 10 to 
60 times smaller than the 7° beam of the COBE DMR. Thus the RAM and CPU time required 
to invert a differential map by this technique would increase by a factor of 10^ to 2 x 10^ over 
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the time required for the COBE DMR maps. Because of this, several groups have abandoned the 
wide-angle differential radiometer design that worked so well in the COBE DMR, and are using 
either total power radiometers or chopping against a load, both of which are prone to serious 
systematic effects. However, we have found a different way of executing the operations in Equation 
^ that reduces the computational load for small beams by a substantial factor. 



4. The Time Ordered Alternative 

Many iterative methods for solving AX = B only require the ability to multiply a vector 
times A. This multiplication can be performed by a sum over the elements of A, or it can be 
performed by scanning the data in time-order and processing each differential pixel pair as it 
occurs. The value of an off-diagonal element in A is just the number of times a given pixel pair 
occurs, and if the typical off-diagonal element is large then using the sparse matrix is much faster 
than scanning the time-ordered data. But the typical value of an off-diagonal element scales like 
Nt/Noff, where Nt is the total number of data samples taken. The proposed new experiments are 
all at least 10 times more sensitive per sample than the COBE DMR, so Nt will be only somewhat 
larger than the COBE DMR data set, while Ngfj will be at least 10^ times larger than it was for 
DMR. Thus for maps with a large number of pixels, the typical value of an off-diagonal element of 
the matrix will be quite small, so the time-ordered approach, which scales like Nt, will be faster 
than scanning the sparse matrix, which scales like N^ff. In addition, the RAM required by the 
time-ordered technique scales like Npi^ while storing the sparse matrix scales like N^jf- Since new 
experiments with large maps will have data rates only somewhat larger than COBE^s, and since 
the speed of computers has grown by orders of magnitude, producing large all-sky maps from 
differential data is easily feasible for future experiments. 

The iterative scheme in Equation ^ can be implemented in time order using the following 
scheme 









- -sf*)^ 
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J2t {^i,p{t) + ^iMt)] 
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The denominator in Equation |6| is just the diagonal of A, and does not need to be computed 
during each iteration. For each data point this algorithm requires obtaining the time-ordered data 
S{t), p{t), and m{t) from disk or tape, two fetches, four adds and two stores. It goes very much 
faster than the actual data-taking. 

In effect, this scheme evaluates the temperature in each map pixel by averaging all the 
differential observations of that pixel, after correcting each for an estimate of the signal in the 
reference beam, obtained from the previous iteration. We have tested this algorithm by simulating 
one year of space mission observations using a DEC 3000/600 Alpha workstation. The input 
map, shown in Figure |2|, contains 1,572,864 pixels and includes both simulated CMB anisotropy, 
including the dipole, and a Galactic foreground model. We "observed" this map with differential 
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Fig. 1. — Raw differential data from simulated space mission observations of the sky shown in 
Figure ||. The dominant signal is the CMB dipole modulated by the motion of the spacecraft. 
Galactic plane crossings are seen as alternating spikes while the small amplitude fluctuations are 
from a cold dark matter CMB simulation. 



radiometers with beam separation 6chop = 135° on a spinning spacecraft. Figure || shows a sample 
of the time-ordered data from this simulation, using a data rate of 20 observations sec~^. We 
processed one year's worth of differential observations, including instrument noise, solving for 
the map using Equation The iterations were started with a pure dipole map for X^^^ . Figure 
^ shows the input map, and the iterations X^^\X^^^ and X^'^^\ Note the large artifacts in the 
first iteration due to the reference beam crossing the Galactic plane, which is not modeled in 
X^^\ These artifacts are quickly removed in subsequent iterations, and by = 20 there are no 
artifacts remaining in the map above the level of the instrument noise. Each iteration requires 
approximately 8 hours of CPU time on the existing Alpha workstation. The required time is 
proportional to the data rate and is virtually independent of the number of map pixels: increasing 
the number of pixels by a factor of 4 slows processing by only 2%. Significantly, 64% of the CPU 
time is spent in the subroutine that simulates the spacecraft attitude. Thus, if applied to real data 



Fig. 2. — Top left: Full sky map used as input for the mission simulation. The map includes 
simulated CMB anisotropy, the CMB dipole, and a model Galactic signal. Top right: The pure 
dipole signal used for the O*'^ iteration. Bottom left: The recovered map after one iteration of 
Equation ^. The Galactic plane signal appears coherently, though with echos that are ~ 10% of 
the plane signal. This efficient reduction of plane echoes after only one iteration requires a scan 
pattern that successfully connects a given sky pixel to many other pixels, both on and off the 
Galactic plane. Bottom right: The recovered map after 20 iterations, by which time no significant 
artifacts remain. 

with stored attitude information, the processing would be even faster. 

This technique can be readily extended to include observations of microwave polarization by 
accounting for the orientation of the observing horns with respect to the sky. We have tested this 
extension with similar mission simulations and find suitable convergence on maps of the Stokes 
parameters I, Q, and U . The computational cost of including polarization is less than a factor of 
2, relative to the pure intensity case. 



5. Discussion 

We have demonstrated a computational method that allows for the production of mega-pixel 
CMB anisotropy maps from large differential radiometer data sets. Thus the proposed new 
anisotropy experiments, which would produce maps with lOO's of times more pixels than COBE 
DMR's, can still use the time-tested wide-separation differencing technique successfully employed 
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by the DMR. 
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Fig. 3. — Same as Figure |2| but with the dipole removed and the temperature scale expanded. The 
O*'^ iteration is omitted from this panel. 



